Computing the monomer-dimer systems through matrix permanent 
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The monomer-dimer model is fundamental in statistical mechanics. However, it is #P-complete 
in computation, even for two dimensional problems. A formulation in matrix permanent for the par- 
tition function of the monomer-dimer model is proposed in this paper, by transforming the number 
of all matchings of a bipartite graph into the number of perfect matchings of an extended bipartite 
graph, which can be given by a matrix permanent. Sequential importance sampling algorithm is 
applied to compute the permanents. For two-dimensional lattice with periodic condition, we ob- 
tain 0.6627 ± 0.0002, where the value is ha = 0.662798972834. For three-dimensional problem, our 
numerical result is 0.7847 ± 0.0014, which agrees with the best known bound. 

PACS numbers: 05.50.+q, 02.10.Ox, 02.70.Uu, 02.50.-r 



I. INTRODUCTION 

The monomer-dimer model is considered, in which the 
set of sites in a lattice is covered by a non-overlapping 
arrangement of monomers (molecules occupying one site) 
and dimers (molecules occupying two sites that are neigh- 
bors in the lattice). It is fundamental in lattice statisti- 
cal mechanics 0, Q . A two dimensional monomer-dimer 
model with size m = (mi, 7712) is a rectangle lattice with 
mi x rri2 sites. The two dimensional monomer-dimer sys- 
tems are used to investigate the properties of adsorbed 
diatomic molecules on a crystal surface [3j; the three- 
dimensional systems occur classically in the theory of 
mixtures of molecules of different sizes H as well as the 
cell cluster theory of the liquid state @- More com- 
plete description of the history and the significance of 
monomer-dimer model can be found in [l[ and the refer- 
ences therein. 

All possible monomer-dimer coverings for a given lat- 
tice defines the configuration space of a monomer-dimer 
model. A fundamental question for such a statistical 
mechanics model is to determine the cardinal number of 
the configuration space. Practically, most of the thermo- 
dynamic properties of physical systems can be obtained 
from the number of all possible ways that a given lattice 
can be covered. Thus a considerable attention has been 
devoted to such a counting problem. For a d-dimensional 
cubic lattice with size m — (mi, 777,2, • • • , m^), this cardi- 
nal number is denoted by Z(m,d). It is proved that the 
following limit exists 



hd = lim 



log Z(m, d) 
mim 2 • • ■ md 



The limit hd is called monomer-dimer constant @. 

Even for the simplest two dimensional models, there 
are very few closed form results on the monomer-dimer 
constant. Baxter and Gaunt gives estimates of the 



constants using the asymptotic expansions. Hammersley 
and Menon [9| estimate the hi by calculating lower and 
upper bounds. Numerical simulation should play a very 
important role. However it has been proved that comput- 
ing the monomer-dimer constant is a #P-complete prob- 
lem even for 2-dimensional problems [loj . which shows 
the hardness of the computation. The Monte Carlo 
method is applied to study the problem in 0, |(| ITlj |. 
which is a natural consideration. Recently, Friedland and 
Peled [IH give a complete up-to-date theory of the com- 
putation of monomer-dimer constant by calculating lower 
and upper bounds. They obtain h 2 = 0.66279897, which 
agrees with the heuristic estimation e h2 = 1.940215351 



due to Baxter 



and 0.7653 < ha < 0.7862. Two- 
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dimensional model with fixed dimer density is studied 
intensively by Kong [13]. The monomer-dimer constant 
with 12 digits accuracy for two dimensional problem is 
given as h 2 = 0.662798972834. 

In this paper, we propose a formulation that trans- 
forms the counting of all matchings of a bipartite graph 
to the counting of perfect matchings of an extended bi- 
partite graph. Hence, the monomer-dimer constants in 
any dimensions can be computed by permanents of matri- 
ces. Permanent of matrix is studied for a very long time 
[hi , Il5j |. After Valiant proves that evaluating the per- 
manent of a 0-1 matrix is a #P-complete problem [16|, 
many randomized approximate algorithms are developed 
[l7l Il8l [l9| . They can give reasonable estimations for 
permanent within a acceptable computer time. 

We consider cubic lattice with periodic condition, and 
concentrate on two and three dimensional lattices in the 
computation. The algorithms are applicable to other di- 
mensions and domains other than rectangle. For simplic- 
ity of notation, we assume that mi = 7772 = • • • = md- 
But this is not essential for the algorithms. 

In the next section, the formulation of the monomer- 
dimer configuration space in matrix permanent is pre- 
sented. Computational methods are discussed in section 
III. The sequential importance sampling algorithms are 
given to compute matrix permanents. In section IV, nu- 
merical results are presented which clearly shows the ef- 
ficiency of our formulation and the computational meth- 
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ods. Finally in section V, some discussions and comments 
are given. 

II. FORMULATION IN PERMANENT 

Consider each point /site in the lattice as a vertex, and 
an edge exists if the two vertices are neighbors in the 
lattice. Hence a graph G — (V, E) is naturally defined. 
Using the terminology of graph theory a monomer-dimer 
system can be represented as a covering of the vertices of 
the graph G = (V, E) by a non-overlapping arrangement 
of monomers (molecules covering one vertex) and dimcrs 
(molecules covering a pair of adjacent vertices). 

It is convenient to identify monomer-dimer configura- 
tions with matching in the graph G. The sites of a cubic 
lattice can be divided into two vertex sets V\ and V 2 . A 
site and its neighbor should always belong to different 
vertex sets. There are edges between neighbors, and all 
edges form a edges set E. Thus an undirected bipartite 
graph G(Vi UV2, E) is constructed. In terms of the graph 
theory, a covering of all vertices with dimers is a per- 
fect matching of the bipartite graph G(V\ U V2,E); and 
a covering with k dimers is a k- matching of it. Hence 
the cardinal number of the configuration space of the 
monomer-dimer model equals to the number of all possi- 
ble matchings of the bipartite graph G. 

The partition function of the system is defined as 

a 

Z(A) = Z G (A)=^m fe A fe (1) 

k=Q 

where rrik = m,k{G) is the number of k— matching in the 
graph G, which is equivalent to the number of monomer- 
dimer configurations with k dimers. Zq(1) enumerates 
all possible matchings in G. 

Let G be a bipartite graph and A be the adjacent ma- 
trix of the graph G. The number of perfect matchings 
of G is equal to the permanent of the matrix A, which is 
defined as 

n 

perm,{A) = ^ JJa^). (2) 

o-en„ i=i 

Here n„ is the symmetric group of degree n. 

A matrix permanent formulation for enumerating k- 
matching of a bipartite graph is proposed by Friedland 
and Levy recently [2l[. Their method can be applied to 
approximate the monomer-dimer constant. For any given 
k, method by Friedland and Levy can compute mk, the 
number of fc-matching, for all k E {0, 1, ■•■ , n}. Thus 
Zq{\) 1 all possible matchings in G, can be given by 

n 

Z G (l) = Y,m k . (3) 

Note that the number of matrix permanents computed 
is n, and n would not be a small number. Here in the 



following we propose a new formulation in matrix per- 
manent. The number of all possible matchings, that is 
Zg(1) in ([3]), can be approximated directly. 

Let A be the adjacent matrix of a bipartite graph G. 
Thus A is a 0-1 matrix. We use G(A) denote the bipartite 
graph with adjacent matrix A. The vertex set of G(A) 
is denoted as V = Vt U V% and the edge set is E. An 
auxiliary graph is constructed based on graph G(A) as 
follows. Vertex sets V 1 and V2 are added to V\ and V2 
respectively. The cardinal numbers of the new sets V 1 
and V-2 are both n. There are n edges between V\ and 
V 2 and each vertex in V\ is adjacent to a different vertex 
in V 2 ■ The vertexes of V 1 are adjacent to every vertexes 
of V2 and V 2 ■ Let 

B=( 1 A f" x "V (4) 

y x nxn x nxn J 

where l nxn is the n x n matrix whose entries are all 
equal to 1; and I n xn is the identity matrix of order n. It 
is obvious that B is a 0-1 matrix, and it is the adjacent 
matrix of the auxiliary graph. 

Let AM (A) denote the number of all possible match- 
ings of the graph G(A). Note that perm(B) gives the 
number of perfect matchings of G(B). In a perfect 
matching of G(B), each vertex V\ is assigned to be ad- 
jacent to a vertex in V% \] V 2 ■ The number of all the 
possible assignment between V\ and V2 \J V 2 equals to 
AM (A). If the adjacent edges between the set V\ and set 
V 2 1J V 2 are chosen, there are nl possibilities for choosing 
the adjacent edges between V x and V2 1J V 2 . So we have 
AM {A) ■ n\ = perm(B), that is 

AM (A) = 1 perm ( A | nxrl \ . (5) 

n\ \ tnxn t nX n J 

Denote 

f{\) = -.verm{ A \ W ) = £ A A*. 

v 7 k=0 

Let mk = rrik(G(A)) be the number of k— matching in 
the graph G(A). It is easy to verify that, 

m k = /„_fc. (6) 

Thus we can get the following permanent formulation 
of the partition function of monomer-dimer system. 

a 

Z{\)=Z G {\)=Y J fn-k\ k . (7) 
k=0 

Hence the partition function of the monomer-dimer sys- 
tem is formulated as matrix permanent. It is important 
to notice that the matrix B is very special in structure, 
which will be explored in the following numerical algo- 
rithms. 
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III. COMPUTATIONAL METHODS THROUGH 
PERMANENT 

Matrix permanent is a long-studied mathematical 
problem in its own right [3, EH- A bridge between 
the computation of permanent and monomer-dimer con- 
stant is established via the relationship ([5]). Thus the 
monomer-dimer constant can be computed by taking the 
advantage of the efficient algorithms in matrix perma- 
nent. 

The definition of the permanent perm(A) looks sim- 
ilar to that of the determinant det(A). However it is 
much harder to be computed. Valiant [la ] proves that 
computing a permanent is a #P-complete problem in 
counting, even for 0-1 matrices. Hence approximate al- 
gorithms, which can give a reasonable estimation for 
perm(A) within acceptable computer time, attract much 
attentions recently. 

Practical approximate methods for matrix permancnts 
are Monte Carlo algorithms. One way to do so is to re- 
late matrix permanents to matrix determinants by ran- 
domizing the elements of matrices [l?], EH ■ The Markov 
chain Monte Carlo approach can give a fully polynomial 
randomized approximation scheme for the permanent of 
any arbitrary nonnegative matrix. This is obtained by 
M.Jerrum, A.Sinclair and E.Vigoda pjj ]. Beichl, O'Leary 
and Sullivan 11| compute the number of /c-matching of 
monomer-dimer model using Markov chain Monte Carlo 
method. They improve the KRS method Q. 

The Monte Carlo methods with sequential importance 
sampling, which are a kind of efficient algorithms for ap- 
proximating permanent, seem to be promising for the 
monomer-dimer problem [2(| H2, HH ■ Beichl and Sulli- 
van give the best known numerical result for 3-D dimcr 
constant by using the techniques 23] . The framework of 
sequential importance sampling for the permanent of a 
0-1 matrix A is as follows. 

Algorithm SIS-P 

Step 1. Choose a nonzero element from the first row 
of the matrix A with some probability p\ . Suppose the 
column index of this element be k\. Set all the other 
entries in the first row and the fcith column to 0's; 
Step 2. Proceed to the next row, applying the same 
sampling strategy as step 1 recursively. Hence the values 
Vii • ■ • ,Pn can be obtained; 
Step 3. Compute X 



j_ 

in 



J_ 

pi 



The output X of Algorithm SIS-P is a random vari- 
able. It is an unbiased estimator to the permanent of 
0-1 matrix A. Different strategies of choosing the proba- 
bility distributions would lead to different sequential im- 
portance sampling algorithms. 

Now we apply Algorithm SIS-P to compute the per- 
manent of the matrix B in (|4|). The matrix structure 
is so special that all the elements in the (n + l)th to 
(2n)th rows of B are 1. Hence only the first n rows of 
B are needed to be considered. Assume that one sam- 



TABLE I: Comparison of three SIS algorithms for small 2- 
dimensional lattice, m denotes the size of (m, m) lattice. 
Every algorithm sample 10,000 samples. Value denotes the 
approximate cardinal number of configuration space of the 
lattice, and computer times are given in seconds. 
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41280 4.14 


41225 5.55 


41985 19.47 


41025 



FIG. 1: The lattice is a (4,4) lattice and thus the adjacent 
matrix is 16 x 16. The x-axis denotes the number of samplings 
and the y-axis denotes the error of the approximate cardinal 
number of configuration space. 
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If N samples are obtained by Algorithm SIS-P, the num- 
ber of all matchings can be approximated by 



AM {A) = 



perm(B) 



hi':' '/V ''"$ Y 



Three different importance sampling methods Ras by 
[IH, Liu by [13], and B-S by [H[ are used respectively 
to compute the number of the cardinal number of the 
configuration space. The results are given in Table [U 
The convergence rates of the three algorithms for m = 4 
are also shown in FIG. [TJ Simple examples show that 
both Liu and B-S give good results, and Liu runs faster 
in the computation of monomer-dimer constant. 

According to the law of large numbers, the mean value 
of these samples gives an approximation to the perma- 
nent. But in fact, the number of samples in our com- 
putation is not really "large" . More precisely, a typical 
sample number in our computation would be 100, 000, 
while the cardinal number of the sample space could be, 
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for example, 10 115 (the two dimensional monomer-dimer 
model with m = 20). 

Notice that the probability distribution of the random 
variable Y = log A looks similar to the normal distribu- 
tion. If the probability distribution of Y is normal with 
A(/i, c), then the expectation of X would be 

E(X) = E(e Y ) = e"+^. 

Other than computing the sample mean of X directly, 
we can estimate the sample mean fl and sample standard 
deviation a of the random variable Y first. 

IV. EXPERIMENTAL RESULTS FOR 
PERIODIC LATTICES 



that the running times for both SIS and A-PRE grow 
polynomially with respect to m. The time complexity of 
SIS, the method developed in this paper, is about 0(m 3 ) 
for 2-dimensional lattice, while the A-PRE, the MCMC 
method by [ill ] . is about 0(m 5 ). Hence it is easy to tell 
that the algorithm SIS is faster distinctly. This suggests 
that the method SIS can be applied to large monomer- 
dimer problems. 

FIG. 2: The relations between the running time of A-PRE 
method and SIS method with the lattice size m are shown. 
The times of SIS method are costs of 100, 000 samples, those 
of A-PRE method are taken from [ill ]. 



The algorithm SIS is used to approximate permanents, 
which gives approximation to the monomer-dimer con- 
stants. The algorithms are programmed in Matlab 7.0 
and all computations in this paper run on Dell PC with 
CPU 2.8G Hz. 



A. Experiments on two dimensional lattices 

Computational results for 2-dimensional monomer- 
dimer problems with periodic boundary conditions are 
presented in TABLE [H] 



TABLE II: m denotes the size of the planar (m, m) lattice. 
Every time, we sample 100, 000 samples and compute the ap- 
proximate result of logZ(m, 2)/m 2 . We do this several times. 
SIS gives the median value of the approximate values; Time 
denotes the time in second for one sampling. A-PRE is the 
value given in [ill ] 
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In order to fit the limit of logZ(m,2)/m 2 as m goes 
to infinity, we apply regression to the computed mean 
values. The regression function is the same as [23| 



y = — 



P2, 



(8) 



where x denotes the lattice size m, y denotes the h2(m) 
and P2 is the monomer-dimer constant. The monomer- 
dimer constant of 2-dimcnsional problem with periodic 
boundary can be obtained from the regression 

h 2 = 0.6627 ± 0.0002 with 95% confidence. 

The approximate results of the monomer-dimer constant 
coincides with the value h 2 = 0.662798972834 by [H 
very well. 



Let compare with the results and the algorithm A- 
PRE, a Markov Chain Monte Carlo method used by 
Beichl, O'Leary, and Sullivan [111]. Though computers 
used here are different, one can still tell the trends in 
the running times. The curve fitting results for algo- 
rithms SIS and A-PRE are shown in FIG [2] It is clear 



B. Experiments on three dimensional lattices 

For 3-dimensional problem with periodic condition, 
computational results are shown in TABLE IThI The time 
complexity for algorithm SIS in 3 dimensional problems 
is about 0(m 6 ). 
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TABLE III: m denotes the size of the cubic (m, m, m) lattice. 
Every time, we sample 100, 000 samples and compute the ap- 
proximate result of logZ(m, 3)/m 3 . We do this several times. 
SIS gives the median value of the approximate values; Time 



denotes the time in 


second for one sampling 


;. A-PRE is the 


value given in [ill] 






m SIS 


Time(sec) 


A-PRE 


4 0.787359 


0.0039 


0.7844 


6 0.786661 


0.0082 


0.7847 


8 0.785821 


0.0345 


0.7870 


10 0.787093 


0.0919 




12 0.785054 


0.2483 




14 0.783476 


0.6693 





To fit the limit of logZ{m, 3)/to 3 as m goes to infinity, 
we apply regression again. The function we use is 

y=-+P2, (9) 

x 

where x denotes the lattice size m, y denotes the h 3 [ra) 
and p2 is the monomer-dimer constant. The result is 

h 3 = 0.7847 ± 0.0014 with 95% confidence. 

This agrees well with the best known bound 0.7653 < 
h 3 < 0.7862 



V. DISCUSSIONS AND COMMENTS 

The construction of the auxiliary bipartite graph is the 
key step in our formulation. Hence the permanent of the 
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